home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libblas / src_original / zgbmv.f < prev    next >
Encoding:
Text File  |  1994-08-02  |  12.0 KB  |  386 lines
  1. *
  2. ************************************************************************
  3. *
  4. *     File of the COMPLEX*16       Level-2 BLAS.
  5. *     ==========================================
  6. *
  7. *     SUBROUTINE ZGEMV ( TRANS, M, N, ALPHA, A, LDA, X, INCX,
  8. *    $                   BETA, Y, INCY )
  9. *
  10. *     SUBROUTINE ZGBMV ( TRANS, M, N, KL, KU, ALPHA, A, LDA, X, INCX,
  11. *    $                   BETA, Y, INCY )
  12. *
  13. *     SUBROUTINE ZHEMV ( UPLO, N, ALPHA, A, LDA, X, INCX,
  14. *    $                   BETA, Y, INCY )
  15. *
  16. *     SUBROUTINE ZHBMV ( UPLO, N, K, ALPHA, A, LDA, X, INCX,
  17. *    $                   BETA, Y, INCY )
  18. *
  19. *     SUBROUTINE ZHPMV ( UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY )
  20. *
  21. *     SUBROUTINE ZTRMV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX )
  22. *
  23. *     SUBROUTINE ZTBMV ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX )
  24. *
  25. *     SUBROUTINE ZTPMV ( UPLO, TRANS, DIAG, N, AP, X, INCX )
  26. *
  27. *     SUBROUTINE ZTRSV ( UPLO, TRANS, DIAG, N, A, LDA, X, INCX )
  28. *
  29. *     SUBROUTINE ZTBSV ( UPLO, TRANS, DIAG, N, K, A, LDA, X, INCX )
  30. *
  31. *     SUBROUTINE ZTPSV ( UPLO, TRANS, DIAG, N, AP, X, INCX )
  32. *
  33. *     SUBROUTINE ZGERU ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA )
  34. *
  35. *     SUBROUTINE ZGERC ( M, N, ALPHA, X, INCX, Y, INCY, A, LDA )
  36. *
  37. *     SUBROUTINE ZHER  ( UPLO, N, ALPHA, X, INCX, A, LDA )
  38. *
  39. *     SUBROUTINE ZHPR  ( UPLO, N, ALPHA, X, INCX, AP )
  40. *
  41. *     SUBROUTINE ZHER2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA )
  42. *
  43. *     SUBROUTINE ZHPR2 ( UPLO, N, ALPHA, X, INCX, Y, INCY, AP )
  44. *
  45. *     See:
  46. *
  47. *        Dongarra J. J., Du Croz J. J., Hammarling S.  and Hanson R. J..
  48. *        An  extended  set of Fortran  Basic Linear Algebra Subprograms.
  49. *
  50. *        Technical  Memoranda  Nos. 41 (revision 3) and 81,  Mathematics
  51. *        and  Computer Science  Division,  Argonne  National Laboratory,
  52. *        9700 South Cass Avenue, Argonne, Illinois 60439, US.
  53. *
  54. *        Or
  55. *
  56. *        NAG  Technical Reports TR3/87 and TR4/87,  Numerical Algorithms
  57. *        Group  Ltd.,  NAG  Central  Office,  256  Banbury  Road, Oxford
  58. *        OX2 7DE, UK,  and  Numerical Algorithms Group Inc.,  1101  31st
  59. *        Street,  Suite 100,  Downers Grove,  Illinois 60515-1263,  USA.
  60. *
  61. *
  62. ************************************************************************
  63. *
  64.       SUBROUTINE ZGBMV ( TRANS, M, N, KL, KU, ALPHA, A, LDA, X, INCX,
  65.      $                   BETA, Y, INCY )
  66. *     .. Scalar Arguments ..
  67.       COMPLEX*16         ALPHA, BETA
  68.       INTEGER            INCX, INCY, KL, KU, LDA, M, N
  69.       CHARACTER*1        TRANS
  70. *     .. Array Arguments ..
  71.       COMPLEX*16         A( LDA, * ), X( * ), Y( * )
  72. *     ..
  73. *
  74. *  Purpose
  75. *  =======
  76. *
  77. *  ZGBMV  performs one of the matrix-vector operations
  78. *
  79. *     y := alpha*A*x + beta*y,   or   y := alpha*A'*x + beta*y,   or
  80. *
  81. *     y := alpha*conjg( A' )*x + beta*y,
  82. *
  83. *  where alpha and beta are scalars, x and y are vectors and A is an
  84. *  m by n band matrix, with kl sub-diagonals and ku super-diagonals.
  85. *
  86. *  Parameters
  87. *  ==========
  88. *
  89. *  TRANS  - CHARACTER*1.
  90. *           On entry, TRANS specifies the operation to be performed as
  91. *           follows:
  92. *
  93. *              TRANS = 'N' or 'n'   y := alpha*A*x + beta*y.
  94. *
  95. *              TRANS = 'T' or 't'   y := alpha*A'*x + beta*y.
  96. *
  97. *              TRANS = 'C' or 'c'   y := alpha*conjg( A' )*x + beta*y.
  98. *
  99. *           Unchanged on exit.
  100. *
  101. *  M      - INTEGER.
  102. *           On entry, M specifies the number of rows of the matrix A.
  103. *           M must be at least zero.
  104. *           Unchanged on exit.
  105. *
  106. *  N      - INTEGER.
  107. *           On entry, N specifies the number of columns of the matrix A.
  108. *           N must be at least zero.
  109. *           Unchanged on exit.
  110. *
  111. *  KL     - INTEGER.
  112. *           On entry, KL specifies the number of sub-diagonals of the
  113. *           matrix A. KL must satisfy  0 .le. KL.
  114. *           Unchanged on exit.
  115. *
  116. *  KU     - INTEGER.
  117. *           On entry, KU specifies the number of super-diagonals of the
  118. *           matrix A. KU must satisfy  0 .le. KU.
  119. *           Unchanged on exit.
  120. *
  121. *  ALPHA  - COMPLEX*16      .
  122. *           On entry, ALPHA specifies the scalar alpha.
  123. *           Unchanged on exit.
  124. *
  125. *  A      - COMPLEX*16       array of DIMENSION ( LDA, n ).
  126. *           Before entry, the leading ( kl + ku + 1 ) by n part of the
  127. *           array A must contain the matrix of coefficients, supplied
  128. *           column by column, with the leading diagonal of the matrix in
  129. *           row ( ku + 1 ) of the array, the first super-diagonal
  130. *           starting at position 2 in row ku, the first sub-diagonal
  131. *           starting at position 1 in row ( ku + 2 ), and so on.
  132. *           Elements in the array A that do not correspond to elements
  133. *           in the band matrix (such as the top left ku by ku triangle)
  134. *           are not referenced.
  135. *           The following program segment will transfer a band matrix
  136. *           from conventional full matrix storage to band storage:
  137. *
  138. *                 DO 20, J = 1, N
  139. *                    K = KU + 1 - J
  140. *                    DO 10, I = MAX( 1, J - KU ), MIN( M, J + KL )
  141. *                       A( K + I, J ) = matrix( I, J )
  142. *              10    CONTINUE
  143. *              20 CONTINUE
  144. *
  145. *           Unchanged on exit.
  146. *
  147. *  LDA    - INTEGER.
  148. *           On entry, LDA specifies the first dimension of A as declared
  149. *           in the calling (sub) program. LDA must be at least
  150. *           ( kl + ku + 1 ).
  151. *           Unchanged on exit.
  152. *
  153. *  X      - COMPLEX*16       array of DIMENSION at least
  154. *           ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
  155. *           and at least
  156. *           ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
  157. *           Before entry, the incremented array X must contain the
  158. *           vector x.
  159. *           Unchanged on exit.
  160. *
  161. *  INCX   - INTEGER.
  162. *           On entry, INCX specifies the increment for the elements of
  163. *           X. INCX must not be zero.
  164. *           Unchanged on exit.
  165. *
  166. *  BETA   - COMPLEX*16      .
  167. *           On entry, BETA specifies the scalar beta. When BETA is
  168. *           supplied as zero then Y need not be set on input.
  169. *           Unchanged on exit.
  170. *
  171. *  Y      - COMPLEX*16       array of DIMENSION at least
  172. *           ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
  173. *           and at least
  174. *           ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
  175. *           Before entry, the incremented array Y must contain the
  176. *           vector y. On exit, Y is overwritten by the updated vector y.
  177. *
  178. *
  179. *  INCY   - INTEGER.
  180. *           On entry, INCY specifies the increment for the elements of
  181. *           Y. INCY must not be zero.
  182. *           Unchanged on exit.
  183. *
  184. *
  185. *  Level 2 Blas routine.
  186. *
  187. *  -- Written on 22-October-1986.
  188. *     Jack Dongarra, Argonne National Lab.
  189. *     Jeremy Du Croz, Nag Central Office.
  190. *     Sven Hammarling, Nag Central Office.
  191. *     Richard Hanson, Sandia National Labs.
  192. *
  193. *
  194. *     .. Parameters ..
  195.       COMPLEX*16         ONE
  196.       PARAMETER        ( ONE  = ( 1.0D+0, 0.0D+0 ) )
  197.       COMPLEX*16         ZERO
  198.       PARAMETER        ( ZERO = ( 0.0D+0, 0.0D+0 ) )
  199. *     .. Local Scalars ..
  200.       COMPLEX*16         TEMP
  201.       INTEGER            I, INFO, IX, IY, J, JX, JY, K, KUP1, KX, KY,
  202.      $                   LENX, LENY
  203.       LOGICAL            NOCONJ
  204. *     .. External Functions ..
  205.       LOGICAL            LSAME
  206.       EXTERNAL           LSAME
  207. *     .. External Subroutines ..
  208.       EXTERNAL           XERBLA
  209. *     .. Intrinsic Functions ..
  210.       INTRINSIC          DCONJG, MAX, MIN
  211. *     ..
  212. *     .. Executable Statements ..
  213. *
  214. *     Test the input parameters.
  215. *
  216.       INFO = 0
  217.       IF     ( .NOT.LSAME( TRANS, 'N' ).AND.
  218.      $         .NOT.LSAME( TRANS, 'T' ).AND.
  219.      $         .NOT.LSAME( TRANS, 'C' )      )THEN
  220.          INFO = 1
  221.       ELSE IF( M.LT.0 )THEN
  222.          INFO = 2
  223.       ELSE IF( N.LT.0 )THEN
  224.          INFO = 3
  225.       ELSE IF( KL.LT.0 )THEN
  226.          INFO = 4
  227.       ELSE IF( KU.LT.0 )THEN
  228.          INFO = 5
  229.       ELSE IF( LDA.LT.( KL + KU + 1 ) )THEN
  230.          INFO = 8
  231.       ELSE IF( INCX.EQ.0 )THEN
  232.          INFO = 10
  233.       ELSE IF( INCY.EQ.0 )THEN
  234.          INFO = 13
  235.       END IF
  236.       IF( INFO.NE.0 )THEN
  237.          CALL XERBLA( 'ZGBMV ', INFO )
  238.          RETURN
  239.       END IF
  240. *
  241. *     Quick return if possible.
  242. *
  243.       IF( ( M.EQ.0 ).OR.( N.EQ.0 ).OR.
  244.      $    ( ( ALPHA.EQ.ZERO ).AND.( BETA.EQ.ONE ) ) )
  245.      $   RETURN
  246. *
  247.       NOCONJ = LSAME( TRANS, 'T' )
  248. *
  249. *     Set  LENX  and  LENY, the lengths of the vectors x and y, and set
  250. *     up the start points in  X  and  Y.
  251. *
  252.       IF( LSAME( TRANS, 'N' ) )THEN
  253.          LENX = N
  254.          LENY = M
  255.       ELSE
  256.          LENX = M
  257.          LENY = N
  258.       END IF
  259.       IF( INCX.GT.0 )THEN
  260.          KX = 1
  261.       ELSE
  262.          KX = 1 - ( LENX - 1 )*INCX
  263.       END IF
  264.       IF( INCY.GT.0 )THEN
  265.          KY = 1
  266.       ELSE
  267.          KY = 1 - ( LENY - 1 )*INCY
  268.       END IF
  269. *
  270. *     Start the operations. In this version the elements of A are
  271. *     accessed sequentially with one pass through the band part of A.
  272. *
  273. *     First form  y := beta*y.
  274. *
  275.       IF( BETA.NE.ONE )THEN
  276.          IF( INCY.EQ.1 )THEN
  277.             IF( BETA.EQ.ZERO )THEN
  278.                DO 10, I = 1, LENY
  279.                   Y( I ) = ZERO
  280.    10          CONTINUE
  281.             ELSE
  282.                DO 20, I = 1, LENY
  283.                   Y( I ) = BETA*Y( I )
  284.    20          CONTINUE
  285.             END IF
  286.          ELSE
  287.             IY = KY
  288.             IF( BETA.EQ.ZERO )THEN
  289.                DO 30, I = 1, LENY
  290.                   Y( IY ) = ZERO
  291.                   IY      = IY   + INCY
  292.    30          CONTINUE
  293.             ELSE
  294.                DO 40, I = 1, LENY
  295.                   Y( IY ) = BETA*Y( IY )
  296.                   IY      = IY           + INCY
  297.    40          CONTINUE
  298.             END IF
  299.          END IF
  300.       END IF
  301.       IF( ALPHA.EQ.ZERO )
  302.      $   RETURN
  303.       KUP1 = KU + 1
  304.       IF( LSAME( TRANS, 'N' ) )THEN
  305. *
  306. *        Form  y := alpha*A*x + y.
  307. *
  308.          JX = KX
  309.          IF( INCY.EQ.1 )THEN
  310.             DO 60, J = 1, N
  311.                IF( X( JX ).NE.ZERO )THEN
  312.                   TEMP = ALPHA*X( JX )
  313.                   K    = KUP1 - J
  314.                   DO 50, I = MAX( 1, J - KU ), MIN( M, J + KL )
  315.                      Y( I ) = Y( I ) + TEMP*A( K + I, J )
  316.    50             CONTINUE
  317.                END IF
  318.                JX = JX + INCX
  319.    60       CONTINUE
  320.          ELSE
  321.             DO 80, J = 1, N
  322.                IF( X( JX ).NE.ZERO )THEN
  323.                   TEMP = ALPHA*X( JX )
  324.                   IY   = KY
  325.                   K    = KUP1 - J
  326.                   DO 70, I = MAX( 1, J - KU ), MIN( M, J + KL )
  327.                      Y( IY ) = Y( IY ) + TEMP*A( K + I, J )
  328.                      IY      = IY      + INCY
  329.    70             CONTINUE
  330.                END IF
  331.                JX = JX + INCX
  332.                IF( J.GT.KU )
  333.      $            KY = KY + INCY
  334.    80       CONTINUE
  335.          END IF
  336.       ELSE
  337. *
  338. *        Form  y := alpha*A'*x + y  or  y := alpha*conjg( A' )*x + y.
  339. *
  340.          JY = KY
  341.          IF( INCX.EQ.1 )THEN
  342.             DO 110, J = 1, N
  343.                TEMP = ZERO
  344.                K    = KUP1 - J
  345.                IF( NOCONJ )THEN
  346.                   DO 90, I = MAX( 1, J - KU ), MIN( M, J + KL )
  347.                      TEMP = TEMP + A( K + I, J )*X( I )
  348.    90             CONTINUE
  349.                ELSE
  350.                   DO 100, I = MAX( 1, J - KU ), MIN( M, J + KL )
  351.                      TEMP = TEMP + DCONJG( A( K + I, J ) )*X( I )
  352.   100             CONTINUE
  353.                END IF
  354.                Y( JY ) = Y( JY ) + ALPHA*TEMP
  355.                JY      = JY      + INCY
  356.   110       CONTINUE
  357.          ELSE
  358.             DO 140, J = 1, N
  359.                TEMP = ZERO
  360.                IX   = KX
  361.                K    = KUP1 - J
  362.                IF( NOCONJ )THEN
  363.                   DO 120, I = MAX( 1, J - KU ), MIN( M, J + KL )
  364.                      TEMP = TEMP + A( K + I, J )*X( IX )
  365.                      IX   = IX   + INCX
  366.   120             CONTINUE
  367.                ELSE
  368.                   DO 130, I = MAX( 1, J - KU ), MIN( M, J + KL )
  369.                      TEMP = TEMP + DCONJG( A( K + I, J ) )*X( IX )
  370.                      IX   = IX   + INCX
  371.   130             CONTINUE
  372.                END IF
  373.                Y( JY ) = Y( JY ) + ALPHA*TEMP
  374.                JY      = JY      + INCY
  375.                IF( J.GT.KU )
  376.      $            KX = KX + INCX
  377.   140       CONTINUE
  378.          END IF
  379.       END IF
  380. *
  381.       RETURN
  382. *
  383. *     End of ZGBMV .
  384. *
  385.       END
  386.